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We study the generation of high harmonic radiation by Bloch electrons in a model transparent solid 
driven by a strong mid-infrared laser field. We solve the single-electron time-dependent Schrodinger 
equation (TDSE) using a velocity-gauge method [New J. Phys. 15, 013006 (2013)] that is numerically 
stable as the laser intensity and number of energy bands are increased. The resulting harmonic 
spectrum exhibits a primary plateau due to the coupling of the valence band to the first conduction 
band, with a cutoff energy that scales linearly with field strength and laser wavelength. We also 
find a weaker second plateau due to coupling to higher-lying conduction bands, with a cutoff that 
is also approximately linear in the held strength. To facilitate the analysis of the time-frequency 
characteristics of the emitted harmonics, we also solve the TDSE in a time-dependent basis set, 
the Houston states [Phys. Rev. B 33, 5494 (1986)], which allows us to separate inter-band and 
intra-band contributions to the time-dependent current. We find that the inter-band and intra¬ 
band contributions display very different time-frequency characteristics. We show that solutions in 
these two bases are equivalent under an unitary transformation but that, unlike the velocity gauge 
method, the Houston state treatment is numerically unstable when more than a few low lying energy 
bands are used. 

PACS numbers: 42.65.Ky, 42.65.Re, 72.20.Ht 


I. INTRODUCTION 

Since high harmonic generation (HHG) in inert gases 
was first discovered in 1987 mm, it has become one of 
the major research areas in ultrafast atomic physics. In 
three decades of development, HHG has pushed the tech¬ 
nology for creating tunable extreme ultraviolet (XUV) 
pulses into the attosecond regime mm, and has been 
widely used to probe the ultrafast dynamics of atomic 
and molecular, and solid systems mm- Since the HHG 
process is highly nonlinear, the intensity of the generated 
harmonics is typically orders of magnitude lower than the 
driving laser intensity. This limits the number of photons 
per pulse that can be obtained, meaning that applications 
such as pump-probe spectroscopy and lithography using 
HHG are not presently practical using gas-phase sources. 

Recently, Ghimire et al. discovered that high order 
harmonics can also be generated from a bulk crystal m , 
which has opened up the possibility of studying attosec¬ 
ond electron dynamics in materials. Because of the use 
of a high-density target, solid-state HHG has a poten¬ 
tial for high efficiency. In addition, it may be possible to 
engineer the structure of the solid target on a microme¬ 
ter scale, and thereby design periodic structures that en¬ 
hance the macroscopic phase matching H3HH], further 
boosting the number of photons generated. From a fun¬ 
damental point of view, solid-state HHG is also interest¬ 
ing as a potential tool for addressing and understanding 


‘Electronic address: mwu3@lsu.edu 
' Electronic address: gaarde@phys.lsu.edu 


the ultrafast dynamics of electrons in periodic structures. 

The electron dynamics in a solid interacting with an 
electromagnetic field are generally considered to have a 
contribution from both intra-band and inter-band dy¬ 
namical processes [T6Ul8| . as illustrated in Fig. [lj The 
intra-band dynamics involves k -space motion of an elec¬ 
tron along one (or several) specific bands, while the inter¬ 
band dynamics involves electron transitions between dif¬ 
ferent bands. Although this picture has been widely 
adopted in studying dynamics in solids, the mechanism 
for solid-state HHG has still not been well characterized 
in these terms. Ghimire et al. Hsus] proposed that 
laser-driven Bloch oscillations of an electron wave packet 
on a single conduction band could be the source of the 
non-linear current responsible for HHG. This model was 
supported by recent results on THz HHG by Schubert et 
al. |2flj . In contrast, calculations by Vampa et al. [18] 
using a two-band model indicated that the HHG spec¬ 
trum is dominated by the inter-band current and further¬ 
more that the cutoff is limited by the largest band gap. 
Hawkins et al. [21] suggested that higher-lying bands 
should be included in the laser-solid description in order 
to accurately capture the laser-driven electron dynamics. 

A theoretical framework for the interaction of lasers 
with solids has been constructed using both many- 
electron models pIB] EDI 122] and single-electron models 
[Hum [MU!- The many-electron models are based on 
second quantization, together with a reduction in cor¬ 
relation via a Hartree-Fock decoupling scheme, which 
leads to the well-known semiconductor Bloch-Equations 
(SBE) [28]. In the simplest case, the SBE describe 
an ensemble of correlated two-level systems [29]. The 
SBE approach has been used extensively in describing 
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FIG. 1: (Color online) The band structure used in our calcu¬ 
lation and the scheme of the inter-band and intra-band dy¬ 
namics. The intra-band dynamics involves the motion of the 
electron on the same band, while inter-band dynamics de¬ 
scribes the transitions of the electron between different bands. 
We regard the second band as the valence band and the third 
band as the conduction band. 


optical properties of solids, where it successfully de¬ 
scribes many semiconductor optical experiments such as 
pump-probe, four-wave-mixing, and photon echoes [29K 
[35]. The single-electron models, on the other hand, 
treat the solid as a single electron in an effective peri¬ 
odic potential. The laser-solid interaction is then de¬ 
scribed by the laser-driven single-electron motion in this 
effective periodic potential [33]. This single-electron ap¬ 
proach has been very successful in addressing electron 
dynamics in periodic structures such as Bloch oscilla¬ 
tions, Zener tunneling, and Wannier-Stark localization 
in semiconductor superlattices [33, ^55], optical lattices 
[Ml [37] as well as waveguide arrays [35] [39]. We note that 
the single-electron models are conceptually the closest to 
the hugely successful single-active-electron treatment of 
HHG in atomic and molecular gases, which has yielded 
a number of insights into both the mechanism and con¬ 
trol over the harmonic and attosecond pulse generation 
processes mmmm- 

In this paper, we will follow the analogy to optical 
lattices and HHG in gases and work in a single-electron 
framework. We solve the time-dependent Schrodinger 
equation (TDSE) for an electron interacting with a peri¬ 
odic potential using two different numerical approaches, 
which allows us a number of insights into the HHG pro¬ 
cess. In the first approach, we follow the velocity gauge 
treatment of Korbman et al. }4S| in which the wave func¬ 
tion is expanded in a basis of Bloch states, which means 
that many bands are included in the calculation. We 
find that this method is numerically stable with respect 
to increasing both the laser intensity and the number of 
bands considered, but that it does not allow us to sep¬ 
arately consider the intra-band and inter-band electron 
dynamics. Our second approach is to solve the TDSE 
in a time-dependent basis set [4B], the so-called Houston 
states, in which the intra- and inter-band contributions 
can be naturally separated. We show that while the two 
methods are equivalent under a unitary transformation, 
the Houston state treatment becomes numerically unsta¬ 


ble as the number of bands is increased. 

We find that the resulting harmonic spectra exhibit 
both a primary and a secondary plateau, each with a 
cutoff energy that depends linearly on the laser elec¬ 
tric field strength. The primary plateau is dominated 
by inter-band transitions between the valence band and 
the first conduction band. The secondary plateau is due 
to transitions between the valence band and the higher- 
lying conductions bands. This plateau is much weaker 
than the primary plateau at low intensity, but increases 
rapidly with the intensity and eventually merges with 
the primary plateau. Using the Houston state approach, 
we also separately analyze the time-frequency character¬ 
istics of the inter- and intra-band contributions to the 
current and find that they exhibit very different char¬ 
acteristics, in the intensity regime where the spectrum 
is dominated by one primary plateau. We propose that 
this difference could potentially be used to experimen¬ 
tally address which mechanism dominates the solid-state 
HHG process. 

This paper is organized as follows: In Section II, we 
present the theoretical framework and the results of solv¬ 
ing the TDSE using the velocity gauge approach, and in 
Section III, we solve the TDSE in the Houston basis US] 
and analyze the results of section II in terms of inter¬ 
band and intra-band dynamics. Section IV presents an 
analysis of the numerical instability of the Houston basis 
treatment, and Section V presents a brief summary. 


II. TDSE IN THE BLOCH STATE BASIS 


We consider a linearly polarized laser field propagating 
through a thin crystal along the optical axis. We describe 
the laser-solid interaction in one dimension, along the 
laser polarization which lies in the crystal plane. We 
follow the velocity gauge treatment in [45], in which the 
TDSE reads 

ih^ t \m) = (h 0 +H int ) i m), (1) 

where Hq is the field-free Hamiltonian and H m t is the in¬ 
teraction Hamiltonian between the laser and the electron 


=£ + y ( i > 

H in t = —A(t)p. 


( 2 ) 

(3) 


A(t) is the vector potential, and is related to the electric 
field by 


A(t) = - f E(t')dt'. (4) 

p is the momentum operator and V ( x ) is the periodic 
lattice potential. We have employed the dipole approx¬ 
imation A(x,t) ~ A(t) because the wavelengths we are 
interested in are much larger than the lattice constant. 
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According to Bloch’s theorem, the eigenstates of the field- 
free Hamiltonian are the Bloch states 


Ho \&nk) — ^n{k) 1 4>nk) > (5) 

where n is the band index and the eigenvalues e n (k) rep¬ 
resent the dispersion relations of the bands. Each Bloch 
state can be written as a product of a plane wave and a 
function periodic in the lattice spacing do: 

(x\ ct> nk ) = e ikx U nk (x), (6) 

where U nk {x) satisfy 

U n k(x + a 0 ) = U nk (x). (7) 

Because the vector potential is independent of x , the lat¬ 
tice momentum k is still a good quantum number m- 
This means the dynamics of the different lattice momen¬ 
tum channels are independent, and the TDSE can be 
solved independently for each k 05] , 

To solve the TDSE for a specific k 0 , we express the 
wave function in Bloch states 

IV’fcoW) = ^ ) Cnkp jt) l^nfeo) > (®) 

n 

and solve for the time-dependent coefficients C nko ( t ) 

ih^C nko = C nko e n (k 0 ) + ^ E Cn'koPZ’ (9) 

n' 

where the p nn i matrix element is the integration of the 
momentum operator over a lattice cell in space 


Pko = (<t>nko\p\<t>n'ko) 

1 r a ° fhd\ 

= -J ( 10 ) 


Usually, the p matrix is dominated by its tri-diagonal 
matrix elements, which means the transitions to higher 
bands are most likely to happen through successive tran¬ 
sitions between intermediate bands. Finally, we calculate 
the time-dependent laser-induced current as the sum of 
of the current in each of the different fc 0 channels j ko [35] 
where : 

jk 0 = [Re [<V’k 0 |p|V’fco>] + eA (t)} ■ (U) 

m 


The laser pulse we use has a cos 4 envelope in its elec¬ 
tric field, with a full width at half maximum (FWHM) 
pulse duration of 3 optical cycles for all wavelengths. 
We have considered laser wavelengths A between 2 /irn 
and 5 /im, and intensities between lxlO 10 W/cm 2 and 
2xl0 12 W/crn 2 . The harmonic spectrum is calculated 
as the modulus square |j(u;)| 2 of the Fourier transform 
of the time-dependent current in Eq. (11). Before the 
Fourier transform, we multiply j(t) by a time-dependent 
window function in order to suppress the coherence be¬ 
tween population in the conduction and valence bands 


which would otherwise last forever (in our model) and 
would dominate the spectrum in the region around the 
band-gap energy. The window function matches the en¬ 
velope of the laser pulse. 

Throughout the paper, we use a periodic potential 
V(x) = —Vo(l + cos(27ra"/ao)), with Vo = 0.37 and lattice 
constant ao = 8, both in atomic units. This Mathieu- 
type potential leads to a band structure that can be ex¬ 
pressed in terms of Mathieu functions [33]. It has been 
used extensively in the optical lattice community [49U5T] . 
The resulting band structure has a minimum band gap 
of 4.2 eV and is shown in Fig. 0 Unless otherwise speci¬ 
fied, we have used 51 Bloch states in our expansion of the 
wave function, which means that 51 bands are included 
in the calculations for each k value. Since the lowest band 
(band 1) is deeply bound and very flat, we use band 2 as 
the initially populated valence band. We have checked 
that transitions involving band 1 play a negligible role in 
the harmonic generation dynamics. The initial popula¬ 
tion is a small superposition (Afco = 7r/20ao) of Bloch 
states near fco = 0 on the valence band, corresponding 
to a wave function which is initially spatially delocalized 
through-out the solid. This implies that the valence band 
is near “frozen” so that only a small distribution of pop¬ 
ulation near k = 0 can be excited to higher bands, and 
is similar to the initial condition proposed in m , where 
only a small part (about 2% in our case) of valence band 
electrons are excited to the conduction band where they 
undergo laser-driven Bloch oscillations. It is also similar 
to the initial condition used in quantum well and optical 
lattice systems when inducing Bloch oscillations [25] 05}- 
1521 . We note that Bloch oscillations (usually thought of 
as electron motion in fc-space) can indeed be captured by 
our formalism in which different k’s are uncoupled from 
each other. In this formalism, Bloch oscillations are more 
easily conceptualized as charge oscillations in real space, 
rather than in k space |53j . 

Fig-0 shows harmonic spectra for our model system 
calculated using a laser wavelength A = 3.2 /jm and in¬ 
tensity 4.5 x 10 11 W/cm 2 . This corresponds to a Bloch 
frequency ug = Ea 0 /H = 0.78 eV, where E is the electric 
field amplitude. Although this intensity is low compared 
to the experiment in da, it is high enough to generate 
rich nonlinear dynamics. The harmonic spectrum ex¬ 
hibits both a perturbative regime (harmonic order < 10), 
a plateau regime (10 — 30) and a cutoff (~ 30), very sim¬ 
ilar to the general structure of the harmonic spectrum 
generated by atoms |40| [41] . As we will show in the 
following section, the plateau is due to inter-band tran¬ 
sitions between the conduction and valence band. This 
agrees with the prediction in m ■ However, in contrast 
to that paper, we find that harmonics can be generated 
with photon energies well above the minimum and max¬ 
imum band gap energies, as shown in Figs. 0 and 0 

The two curves in Fig. 0 represent the full calculation, 
including all 51 bands, and a calculation in which bands 
4 and 5 (the second and third conduction band) have 
been dynamically excluded. This is done by setting the 
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FIG. 2: (Color online) High harmonic spectra (logarithmic in¬ 
tensity scale) of the laser induced current calculated by solv¬ 
ing the TDSE in velocity gauge. The laser wavelength and 
peak intensity are 3.2 jim and 4.5 x 10 11 W/cm 2 . The red 
solid curve shows the spectrum where 51 bands are used in 
the calculation whereas the blue dashed curve uses the same 
condition as the red curve, except bands 4 and 5 are removed 
from the calculation. The black arrow indicates the minimum 
and maximum band gap energies. 


coefficients of C\k 0 and C$k Q to zero in Eq. Q at each 
time step in solving the TDSE. We note that since the 
transition probability from bands 2 and 3 to bands 6 and 
above is very small, removing bands 4 and 5 makes our 
model effectively a two-band model. The full-calculation 
spectrum exhibits a weak second plateau, about 10 or¬ 
ders of magnitude lower than the main plateau, which is 
absent in the reduced-band calculation. The comparison 
between the two curves yields two insights: (i) the main 
part of the plateau and cutoff region is well described by 
the dynamics involving just the valence and the lowest 
conduction bands, indicating that higher bands play a 
negligible role in the harmonic generation in this wave¬ 
length and intensity regime, and (ii) the second plateau 
is due to contributions to the dynamics involving higher- 
lying bands, predominantly due to transitions between 
bands 4 or 5 and the valence band. 

We next investigate the intensity and wavelength de¬ 
pendence of the first and second plateau and cutoff. 
Fig. [3ja) shows the harmonic yield as a function of laser 
electric field strength for A = 3.2 gm. We see that the 
cutoff of the first plateau increases linearly with field 
strength. This is in agreement with the experimental 
finding in EJJ- Fig. [3jb) shows the same intensity scan 
without bands 4 and 5, and the linearity of the first cutoff 
is revealed for a larger range. Going back to Fig. §a), 
we see that as the intensity increases, the second plateau 
rises and merges with the first plateau consistent with 
what was observed in [54j . Subsequently, the cutoff en¬ 
ergy of this new, longer plateau, also exhibits a linear 
dependence on the laser field strength. Fig. [3](c) explores 
the build-up of the second plateau in more detail. We 
show a line-out of the field-strength dependence of har¬ 
monics 19 and 61 (H19 and H61), which are located in 
the first and second plateau respectively. At the highest 
fields, these harmonics are both in the plateau and there¬ 


fore have similar yields and change only slowly with field 
strength. However, at the lowest field strengths, H61 is 
essentially zero (not shown in Fig. [3jc) ) until it starts 
to increase exponentially with intensity, indicating that 
the second plateau is not independent but rather built off 
of the first plateau. This is consistent with our finding 
above, that the population in the higher bands is built in 
a step-like process based on first populating lower bands. 



FIG. 3: (Color online) (a) Harmonic yield as a function of 
laser field strength for A = 3.2 gm. (b) Same as (a) but 
excluding bands 4 and 5. White dashed lines indicate the 
linear dependence of the cutoff energy on field strength, and 
vertical black lines indicate the minimum and maximum band 
gaps between the valence and conduction bands, (c) Field 
strength dependence of H19 and H61 from (a), (d) Same as 
(a) but using A = 4.0 gm. (e) Wavelength dependence of the 
cutoff scaling coefficient with electric field strength. All yields 
are shown in logarithmic scale. 


Our model predicts similar harmonic generation dy¬ 
namics at other wavelengths. Fig. [3jd) shows the field 
strength dependence of the harmonic spectrum at A = 
4.0 gtm, which also exhibits two plateau regions with two 
different, linear, dependences of the cutoff energy on the 
field strength. 

From our numerical results, we can quantify the scal¬ 
ing of the cutoff energy in the following way. We start 
by writing the cutoff energy £ C utoff in units of the Bloch 
frequency, as: 


^cutoff ~ /3ws, (12) 

then the scaling factor for the first and second plateaus 
are j3\ = 14, /3 2 = 38, respectively. Fig. [3](e) shows 
the scaling factor for 6 different wavelengths for the first 
plateau and suggests that the cutoff also scales approxi¬ 
mately linearly with the wavelength. Thus, in this model, 
the first cutoff energy depends linearly on both the elec- 
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trie field strength and the wavelength: 

Cutoff (X A E, (13) 

where the proportionality constant depends on the band 
structure. The linear dependence on held strength con¬ 
trasts with the (A E) 2 scaling of the cutoff in atomic and 
molecular gases [53|, but agrees well with the prediction 
for a strongly driven two-level system BEIG21 which in 
our case would be £ cu to ff = PveEA/nc , where p vc is the 
momentum operator matrix element between the valence 
and conduction bands at k = 0. The two-level formula 
underestimates our numerical values for (3 by about 10%. 
The scaling of the second cutoff with wavelength is more 
difficult to quantify. Both Fig. [3](a) and (d) suggest that 
a third plateau appears at the highest energies, possi¬ 
bly due to the contribution of bands 6 and 7. Apart 
from concerns about applying our simple model to such 
high intensities, we note that there are as of yet no high 
harmonic experiments with the level of sensitivity that 
would be needed to observe such an effect. 

To conclude this section, we briefly comment on our 
choice of initial condition in which a small fc-state wave 
packet yields an initial wave function which is spatially 
delocalized across the entire (ID) crystal. Other recent 
calculations have considered a different initial condition 
in which the valence band is initially fully populated 
ns eu, which in our model would correspond to an 
initial wave function localized at one particular lattice 
site. In a real insulating material, the filled valence band 
means that all the different electronic states of the va¬ 
lence band are occupied by different electrons. The full 
valence band thus only has meaning in the multi-electron 
context. In a single-electron framework the valence band 
can never be filled in the same way, since we only have 
one electron, which corresponds to a much lower dimen¬ 
sional Hilbert space than the multi-electron wave func¬ 
tion. In this sense, in the single-electron framework the 
solid is modeled like a super-atom with a atomic potential 
that is periodic. What we can choose is only the initial 
wave function for this super-atom. For instance, we can 
choose the initial condition of the super-atom to be a 
Bloch state (few fc’s, spatially delocalized) or a Wannier 
state (many fc’s, spatially localized). As demonstrated re¬ 
cently in rare-gas clusters [58], the delocalization of the 
initial wave packet may affect the HHG process, which 
could also be the case in solids. 


III. HOUSTON STATE BASIS 


with a momentum fco in a static field is built from the su¬ 
perposition of a large number of bands all at the same fco- 
Obviously, in this picture there can be no separation of 
the current into intra-band and inter-band contributions. 

In this section we describe an alternative way to calcu¬ 
late the electron dynamics using a time-dependent basis 
set, the Houston states [55]. As demonstrated in Ap¬ 
pendix A, the two solutions are equivalent, since they 
are related by a unitary frame transformation. In the 
Houston basis, however, we can obtain a separation of 
the induced current into intra- and inter-band compo¬ 
nents. This will allow us to separately explore the time- 
frequency characteristics of the two contributions, and 
show that they exhibit very different emission times. 

The Houston states are best thought of as an adiabatic 
basis in which the lattice momentum that would be fco in 
the absence of a field has the time-dependence: 

fc(t) = fc 0 + ^M. (14) 

n 

By construction they are the instantaneous eigenstates 
of the time-dependent Hamiltonian H(t ): 


H(t) \4>nk 0 {t )) = £n(k{t)) \(/>nko(t)) , (15) 

where H (t) is the Hamiltonian in the same single-electron 
Schrodinger equation as above Eq. (jT]) [T6] except for an 
additional term proportional to A 2 : 


iH dt = 


(p + eA) 2 
2m 


V(x) 


m))- (16) 


Including this term in the Schrodinger equation makes 
the form of the Houston states simpler, but it has no 
effects on the current since the wave function only differs 
by an overall time-dependent phase. In this convention, 
the Houston states are related to the Bloch states with 
lattice momentum fc(i) by |46j 

\4>nk 0 {i)) = e- ieA&/n \Km), (i7) 

where x is the position operator. Expanding the time 
dependent wave function with initial lattice momentum 
fco in Houston states 

IV'(l)) =^2a nko (t) \4>nk 0 (t)) , (18) 

n 

we find equations of motion for the coefficients 


In the previous section, the electron dynamics was de¬ 
scribed in a static basis of Bloch states using the velocity 
gauge interaction. In this picture the time dependence of 
the wave function is due solely to the time dependence of 
the Bloch state coefficients C nk (t). Though computation¬ 
ally convenient, this method provides a time-dependent 
current which is hard to understand at an intuitive level. 
For example, the familiar Bloch oscillation of an electron 


ih 


ddnko 

dt 


^ ] j^£n(fc(l))^nn / 


eE(t)X nn '{k{t)) 




. (19) 

where we have made use of Eq. ©• X nn f is the inter¬ 
band transition matrix element defined by 


i r° d 

X nn f (fc) = T— / U* k — U n : k dx 
ia 0 J 0 Ok 


( 20 ) 
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where X nn i = 0 if n = n!. It is calculated numerically 
using procedures in [1601 . Note that the initial wave func¬ 
tion has a lattice momentum of kg, but we could also 
make it a wave packet as we did in the previous section. 

The time-dependent Houston states describe the elec¬ 
tron dynamics in a moving frame in which the lattice 
momentum is prescribed by the vector potential as in 
Eq. (14). Pictured in k space, one can think of an electron 


wave packet oscillating on each energy band, while at the 
same time some of the amplitude transitioning between 
different bands, corresponding to the intra- and inter¬ 
band dynamics, respectively. The motion on each band 
is governed by the time-dependent dispersion e n {k(t)). 
For the intensity used in Fig. [2j the motion of the wave 
packet in k space samples about 2/3 of the first Brillouin 
zone. 


The total current can be calculated from Eq. (11), us¬ 


ing Eq. (181 for the wave function: 


Jtot — He 
m 


^ ] a nkn a n'k 0 {4>nk(t) |p|/’n'fc(t)) 


■ ( 21 ) 


Since now the system is described in the frame that moves 
along with the field, there is no A term in the expression 
for the current, as opposed to that of Eq. (11). 


Fig.[4]compares spectra calculated in the Houston basis 
and the Bloch basis at three different field strengths, cor¬ 
responding to the electron wave packet in fc-space sam¬ 
pling about 2/3, 3/3, and 4/3 of the first Brillouin zone. 
In the Houston basis, we use the three lowest bands 
shown in Fig. [l] while in the Bloch state basis we use 
as many of the 51 bands that were used to calculate 
the band structure as are necessary for numerical con¬ 
vergence. Of these three bands, only the valence and 
conduction bands (2 and 3) contribute meaningfully to 
the dynamics. The reason for using only three bands in 
the Houston basis will be discussed below in connection 
with the numerical properties of the X matrix elements. 
The initial condition used in both cases is a single k state 
(k = 0) in the valence band where the band gap is the 
smallest. Fig. |4])a) and (b) shows that the agreement 
between the Houston and the Bloch calculations is excel¬ 
lent for the part of the spectrum that is dominated by the 
dynamics in the valence and conduction band only. This 
agreement is expected, since the two wave functions are 
related by a unitary transformation (see Appendix A). 
Since the Houston calculation only includes three bands, 
it cannot be expected to reproduce the high-frequency 
part of the spectrum that in the Bloch calculation is due 
to higher-lying bands (approximately above harmonic or¬ 
der 50 in both Fig. B^) and (c)). It is worth noting, 
though, that at the highest intensity the slope of the 
Houston spectrum matches that of the Bloch spectrum. 
We will comment on this in more detail below. 

One advantage of the Houston basis is that the electron 
dynamics naturally separate into an intra- and inter-band 
contribution, and can be studied separately. In Eq. (21), 


the intra-band contribution to the current involves only 




FIG. 4: (Color online) Comparison of harmonic spectra (log¬ 
arithmic scale) from a three-band Houston basis calculation 
and a 51-band Bloch basis calculation, for three different in¬ 
tensities of a 3.2 /im driving field, (a) 4.5 x 10 11 W/cm 2 , (b) 
1.0 x 10 12 W/cm 2 , and (c) 1.8 x 10 12 W/cm 2 . The initial 
condition used in both calculations is a single k = 0 in the 
valence band. The thin line in (a) and (c) is calculated from 
the conduction band only, in the Houston basis, by ignoring 
inter-band transitions (see text). 


Houston states on the same band (n = n’), whereas the 
inter-band contribution involves transitions between dif¬ 
ferent bands (n ^ n'). 

Jintra = / J \ &nk 0 \ (4 > nk(t) \p\ fink^)) (22) 


^ 1 a nk 0 a n’k a {(t>nk(t) |p| fyn'k^t) ) 

n,n 
.n^n 

(23) 

The intra- and inter-band contributions to the current 
are shown in Fig. [sja), for the Houston spectrum shown 
in Fig. Ba). We find that for the range of intensities 
where the three-band Houston model is applicable, the 
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inter-band contribution to the plateau in the spectrum 
is stronger than the intra-band contributions by several 
orders of magnitude in the plateau regime. This is in 
agreement with the prediction in [IS]. In Fig. [5jb) , we 
show the intra-band contribution from the valence and 
conduction band separately. This is done by plotting sep¬ 
arately the terms in the sum in Eq. ( [22] ). We note that 
the intra-band contribution from the valence band would 
be unphysical for a real insulator in which the valence 
band would be filled. In the Bloch state basis, we have 
performed calculations in which we used a full valence 
band as the initial condition. This suppresses the yield 
of the few lowest harmonics but otherwise leads to a har¬ 
monic spectrum that agrees well with those in Fig. |4][a) 
until approximately harmonic 30. We conclude from this 
that the (overestimated) contribution from driven Bloch 
oscillations in the valence band can be ignored for the 
range of frequencies we are interested in. 

It is interesting to note that the single conduction band 
model used in nnum comes naturally from the Hous¬ 
ton model if one starts with the initial population in the 
conduction band and eliminates inter-band transitions. 
When the inter-band transition matrix X vanishes, the 
total current reduces to the intra-band current expres¬ 
sion used in m, where the current originates in the 
electron motion in the conduction band (see proof in Ap¬ 
pendix B). The resulting spectrum is plotted in Fig. [4] 
(labeled ’’conduction only”) , normalized to the response 
at the fundamental. For the lower intensities, the intra¬ 
conduction-band current spectrum is completely different 
than the full (Bloch basis) spectrum and does not exhibit 
a plateau or cutoff. This is in agreement with the finding 
that the harmonic spectrum is dominated by inter-band 
transitions. It also shows that the clear cutoff that can 
be seen in the conduction band spectrum in Fig. §b) 
is in fact due to the time-dependence of the population 
transfer in and out of the conduction band (through the 
time-dependence of the transition matrix X), rather than 
due to the intra-band dynamics itself. At the highest in¬ 
tensity, Fig. Qc), the slope of the intra-conduction-band 
spectrum agrees very well with that of the full calcula¬ 
tion, but is again lacking a cutoff energy. In Section II, 
we interpreted the extended (secondary) plateau in the 
harmonic spectrum as being due to transitions involv¬ 
ing high-lying bands, when describing the dynamics in 
the Bloch state picture. The result in Fig. |4](c) suggests 
that an alternative interpretation is that the extended 
plateau is due to driven conduction-band Bloch oscilla¬ 
tions traversing the entire Brillouin zone, but that in this 
model the cut-off energy cannot be captured without con¬ 
sidering inter-band transitions mm- 

For the remainder of this section, we return to the 
lower intensity case in which we observe a clear distinc¬ 
tion between inter- and intra-band dynamics. In order 
to investigate the electron temporal dynamics using the 
Houston basis wave function, we perform a wavelet trans¬ 
form of the intra- and inter-band currents. The wavelet 
transform is similar to a windowed Fourier transform and 



FIG. 5: (Color online) (a) The intra-band (red, dashed) and 
inter-band (blue, solid) current (logarithmic scale) of a Hous¬ 
ton model. The low order harmonics mainly come from intra¬ 
band current while the harmonics in the plateau are mainly 
come from the inter-band current, (b) The valence band (or¬ 
ange) and conduction band (purple) contributions to the intra 
current (logarithmic scale). 


provides time-frequency information about f the two con¬ 
tributions. In the wavelet transform we use an order 10 
Gabor wavelet to achieve a balanced resolution in both 
the time and frequency domains. Fig. [6] shows the result¬ 
ing time-frequency profiles for the two different contribu¬ 
tions, which clearly exhibit two distinct types of dynam¬ 
ics. We will discuss these separately in the following. 




FIG. 6: (Color online) The wavelet transform (linear scale) 
of the (a) inter-band and (b) intra-band current of a Houston 
model. The harmonics generated from inter-band dynamics 
mainly happens at the peak of the vector potential and there 
is a clear chirp, while the harmonics generated from intra¬ 
band dynamics mainly happens at the peak of the field and 
don’t have clear chirp. The yield is saturated for the lowest 
frequencies in (b). 
















The time structure of the inter-band current in 
Fig. 6ja) exhibits two emission times that are symmet¬ 
rically placed around the peak of the vector potential in 
each half cycle. These two emission time profiles have 
opposite chirps, and are merged at the cutoff frequency. 
In the momentum picture we described above, these two 
emission times arise from the fact that the lattice mo¬ 
mentum k(t ), and thereby the time-dependent band gap 
e(k(t )), traverses all allowed values twice in a half cycle. 
Since the two emission times contribute to the current 
with about the same strength, our model suggests that 
one or the other must be filtered out to obtain a train 
of identical attosecond pulses from the plateau. In con¬ 
trast, harmonics near the cutoff frequency are generated 
near the time when the band gap e(k(t)) is the largest, 
corresponding to the peaks of the vector potential in our 
case. 

The time structure of the intra-band current Fig. [§ b ), 
on the other hand, follows simply from the time- 
dependence of the band structure e(k(t)), i.e., the cur¬ 
vature of the bands. The emission times correspond to 
those times when the band has the largest curvature (the 
largest rate of change in the group velocity), which corre¬ 
sponds to the zeroes of the vector potential, as shown in 
Fig.[(|b). The time-frequency features of the intra-band 
current can also be seen from the cosine expansion of the 
conduction band in m 

oo 

jintra(i) = ^ D S SHI [(2s - 1 )wt\ (24) 

3 — 1 

where D s is related to the band structure and the 
strength of the field. From this expansion, it is clear 
that all the harmonics are generated at the same time 
and there is no chirp in the generated Held. We propose 
that the differences in the two time-frequency character¬ 
istics could be used as an experimental signature of the 
intra- and inter-band dynamics. 

An analogy can be drawn between the Houston picture 
and a strongly driven two-level system to help understand 
the picture of the inter- and intra-band dynamics that 
naturally emerges from it. In the two-level system, the 
adiabatic states are the instantaneous eigenstates of the 
time-dependent Hamiltonian. The dynamics of the sys¬ 
tem can then be separated into an adiabatic part and a 
diabatic part using the adiabatic states as the basis. The 
adiabatic motion describes the evolution of the system 
along the adiabatic states, whereas the diabatic motion 
describes the transition between the adiabatic states. In 
solids, this same separation is achieved in the Houston 
states, which are the instantaneous eigenstates (adiabatic 
states) of the system, as shown in Eq. |l5| ). Since the 
adiabatic states are time-dependent themselves, the adi¬ 
abatic evolution generates nontrivial dynamics by itself 
(the Bloch oscillation, see Appendix B). Similarly, the 
inter-band dynamics can also be understood as diabatic 
transitions between adiabatic states, the same as in a 
two-level system. Note that besides these two descrip¬ 
tions, a third commonly used description for a two-level 


system is the Floquet states, which are the true eigen¬ 
states of the laser-dressed system. The counterpart for 
this description in a solid is the Floquet-Bloch theory and 
it has been well studied in mm- 


IV. NUMERICAL DIFFICULTY OF THE 
HOUSTON BASIS 


As a matter of practice, numerical models are more 
trustworthy if convergence can be achieved with respect 
to all of the parameters in the model. In the present case, 
the number of bands would seem to be such a parameter, 
along with the number of k points and time step size. In 
calculations using the Houston basis, however, we find 
that the numerical results converge poorly if we include 
more than three bands in our model. In this section we 
discuss the numerical difficulty of including more bands 
in the Houston basis, which raises questions about the 
validity of the Bloch oscillation picture for higher bands. 

As we include more bands in the Houston model, the 
gap between neighboring bands becomes very small, and 
the X matrix elements in Eq. (201 increase rapidly. In 
fact, the X matrix grows approximately exponentially 
as a function of the number of bands for our current 
parameters, as shown in Fig. [t]A). Since the X matrix 


— P11 — P22 — P33 — P44 — P55 



FIG. 7: (Color online) The maximum matrix element of the 
X matrix seems to increase rapidly as we include more bands. 
The diagonal part of the P matrix is approaching discontinu¬ 
ous as the bands get higher. The parameters used here is the 
same as in Fig. [2] 

is responsible for the inter-band transitions, in order to 
have the TDSE numerically converged, we must have a 
time step that is small enough to resolve its largest matrix 
element. Apart from resolving the X matrix, the time 
step must also be small enough to resolve the largest gap 
energy in the band structure. As we increase the number 
of bands, the requirement by the X matrix of the time 
step becomes the dominant one. For example, to include 
7 bands, the X matrix requires a million time steps per 
laser cycle, which is much more than that required to 
resolve the highest energy difference in the system. 

Another numerical difficulty in the Houston basis is 
the near discontinuity in the momentum operator matrix 
elements for higher bands. As we go to higher bands, 
the band structure will have sharp turning points at the 
band center and band edge. These turning points result 
in an almost discontinuous behavior of the momentum 
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operator matrix element as shown in Fig. 0b). This 
discontinuous behavior can also be expected from the 
relationship between the momentum operator and the 
derivative of the band structure Eq. (B6). 

These numerical difficulties suggest that separation of 
the inter-band and intra-band dynamics may not be an 
optimal picture for an electron in the higher bands where 
Zener tunneling between the neighboring bands is so 
large that an electron never performs solely intra-band 
motion. Instead, it tunnels through the avoided crossing 
almost like a free electron m ■ The artificial separation 
of a Bloch oscillation motion from the almost free elec¬ 
tron motion complicates the physical picture as well as 
makes the numerical calculation difficult. This separa¬ 
tion is the underlying cause of the problematic behaviors 
of the p and X matrix shown in Fig. [7] This is also a gen¬ 
eral problem for studying strongly driven systems using 
an adiabatic basis. In these systems, the avoided cross¬ 
ings between the adiabatic states are so small that the 
transitions between them blow up. In those conditions, 
the true eigenstates of the system (the Floquet states), 
may be a better basis to work with, though it is difficult 
to apply these methods to broadband driving pulses. 


also suggested that, in general, the cutoff in the har¬ 
monic spectrum is tied to inter-band dynamics through 
the time-dependence of the population in the valence and 
conduction bands. 

Finally, we showed that in the regime where the intra- 
and inter-bancl dynamics can be clearly separated, they 
have very different time-frequency signatures, and we 
proposed that this could be harnessed to experimentally 
characterize the harmonic generation dynamics. 
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V. SUMMARY 


Appendix A: The unitary transformation of 
solutions in the Bloch and Houston basis 


We have studied high harmonic generation in a model 
transparent solid using a ID single-electron model and 
found that this system presents rich nonlinear dynamics. 
In the laser induced current, we see a high harmonic spec¬ 
trum with multiple plateaus. Using the numerically ro¬ 
bust Bloch basis, we have shown that the primary plateau 
is due to transitions between the valence band and the 
lowest conduction band, whereas the secondary plateau 
and more generally higher frequencies in the spectrum 
are due to contributions from higher lying bands. We 
find that the cutoff of the primary plateau scales linear 
with the field strength, in agreement with current ex¬ 
periment m, and we predict that this cutoff also scales 
linearly with the driving wavelength. 

We have also shown that the dynamics of our model 
system can be expressed in either the Bloch basis or the 
Houston basis, and solutions in these two basis are con¬ 
nected through a unitary transform. The Houston basis 
allow for an intuitive separation of intra- and inter-band 
dynamics, and we found that for moderate intensities 
the harmonic radiation is due primarily to inter-band 
dynamics, in agreement with the prediction in (18] . At 
higher intensities, though, this artificial separation be¬ 
comes more problematic as more bands are strongly cou¬ 
pled to each other which manifests itself in the numerical 
calculation becoming unstable. By limiting the dynam¬ 
ics in the Houston basis to just the conduction band, we 
found that for these high intensities, an alternative inter¬ 
pretation of the extended spectral range is provided by 
driven conduction-band Bloch oscillations that traverse 
the entire Brillouin zone. Our Houston basis calculations 


In this section, we show that the solution of the TDSE 
in the Bloch state basis is connected to that in the Hous¬ 
ton basis by a unitary transformation. In the Bloch state 
approach, we solve the TDSE in the velocity gauge, and 
the equation reads 


ih §~ t IV’feoW) = 


p“ T ^_ N eA(t) 


L +v ^ + 


m 


!<(*))> (Al) 


where fco labels the single lattice momentum channel that 
we consider, since different channels are independent. We 
then express the solution in the basis of Bloch states, 
which are the eigenstates of the field-free Hamiltonian 

IV’fcoW) =J2 C n{t) l</W > ( A2 ) 

n 


where c n ’s are the time-dependent energy band ampli¬ 
tudes that this model solves for. 

In the Houston approach, TDSE reads 



(p + eA(t)) 2 
2m 


V(x) 


!<(*)>• (A3) 


The wave function is expressed in the Houston states 

\4>k 0 (t)) = J2 a nk 0 (t) \4>nk 0 (f)) , (A4) 

n 

The Houston states are related to the Bloch states by 

\KkM = e- ie&A{t)/h \<t>n K t)) ■ (A5) 
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We can then match the wave function in these two basis, 
and get the unitary transformation matrix between the 
expansion coefficients. Since \tp^) and IV^) satisfy the 
Schrodinger’s equations that differ by a time-dependent 
A 1 2 3 4 5 term, they must be related by a phase factor 

KW) = e- a/ X(f)) (A6) 

where 

A = —— (A7) 

2m J o 

Substitude Eq. ( |A2| ) and Eq. ( |A4[ ), we finally come to 

a nko (t) = e^5> fe0 (i) (tn H t)\e ie&A(t)/h \K'k*) 
n' 

(AS) 

which is the unitary transformation we are seeking. Note 
that this unitary transformation is very similar to the 
Kramers-Henneberger transformation in the atomic case, 
where a spatial transformation shifts the system into 
the accelerated frame, corresponding to the motion of 
a charged particle in the electric field [BJ. The Kramers- 
Henneberger transformation is a transformation in space 
whereas the transformation here is in momentum space. 


Appendix B: Connection of the single band model to 
the intra-band motion 


the population on each band stays the same as the initial 
condition. The solutions for Eq. (A3) are exactly the 
Houston states apart from a phase 


a nka (t) = a nko (0)e-ifo^(Ht'))dt' (B1) 

!<(*)> = e~ ieA W h J2*n k0 (t) \Km) • (B2) 


Substituting a n k 0 into Eq. (21) the total current is 


Jtot ^ ^ ( |^n&o(0)| i$n k (t) • (B3) 


If we consider the same situation as in mm where 
initially only the lowest conduction band is populated 
and other bands are empty, then the total current reduces 
to 


g 

Jtot = ~ Ipl^cfcp)) • (B4) 

This then reduces to the single conduction band model 
used in mm, where the semi-classical current is de¬ 
rived from a group velocity 


e de c 

Jtot = -ev g = 


(B5) 


k=k(t ) 


since the diagonal elements of the momentum operator 
is related to the band structure by 


In this section we show that if we prevent the inter¬ 
band dynamics, we essentially come back to the single 
conduction band model used in cdi na¬ 
if the inter-band dynamics is not allowed, then the 
inter-band transition matrix X vanishes in Eq. (19), and 


(<Pnk\p\<t>n k ) = n ^ d ^- ( B6 ) 

The proof of the last step can be found in many text¬ 
books. For example, see chapter III in |64j or Appendix 
Ein|653. 
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